#include "StdAfx.h"
#include "Calculator.h"
#include <algorithm>

//////////////////////////////////////////////////////////////////////////
// Constants
static const int kMagnificationBit = 12;
static const int kMagnification = 1 << kMagnificationBit;
static const int kTwiceMagnification = 2 * kMagnification;
static const int kMagnificationMinusOne = kMagnification - 1;

static const int kInitialBufferLength = 1024;

static const float kPoorApproximationPrecision = 2.0f;
static const float kParabolaThreshold = 0.125f;
static const float kBestApproximationToParabolaFactor = 0.38f;

// these constants should correspond to the items in the submenu
// "Curvature calculation"
enum CurvatureCalculationToken
{
	kCurvatureCalculationNone = 0,
	kCurvatureCalculationFast,
	kCurvatureCalculationBest
};


///////////////////////////////////////////////////////////////////////////
Calculator::Calculator(const IBitmapPicture& inBitmap, float inBlurRadius)
	: mSubpixelInfoProvider(inBitmap, inBlurRadius),
	mCurrentDotBooker(&mAnchorDotBooker),
	mBufferKeeper1(kInitialBufferLength, sizeof(Info)),
	mBufferKeeper2(2 * kInitialBufferLength, sizeof(Info*)),
	mInitialDotHalfsizeBig(int(inBlurRadius * kTwiceMagnification)),
	mThresholdFactor(kParabolaThreshold),
	mRegardEtaVariations(false)
{
	SetMaximumDotSize(4 * inBlurRadius);
}


Calculator::~Calculator()
{
}

const SubpixelInfo* Calculator::GetInfoBuffer() const
{
	return mSubpixelInfoProvider.GetInfoBuffer();
}

const SubpixelInfo* Calculator::GetInfo(int inX, int inY) const
{
	return mSubpixelInfoProvider.GetInfo(2 * inX, 2 * inY);
}


void Calculator::SetThresholdFactor(float inThresholdFactor)
{
	mThresholdFactor = inThresholdFactor;
}

void Calculator::SetInitialBounds(float inHalfsize, float inCurvature)
{
	mMaxXi = mMaxEta = mMaxLimit = int(inHalfsize * kTwiceMagnification);
	mMinXi = mMinEta = - mMaxXi;
	mScaledCurvature = inCurvature / kTwiceMagnification;
}

void Calculator::SetInitialDefaultBounds()
{
	mMaxXi = mMaxEta = mMaxLimit = mInitialDotHalfsizeBig;
	mMinXi = mMinEta = - mMaxXi;
	mRegardEtaVariations = false;
	mScaledCurvature = 0.0f;
}

void Calculator::ExtendXiBounds()
{
	int deltaXi = int(mDeltaXi);
	if (deltaXi > 0)
	{
		mMaxXi += deltaXi;
		if (mMaxXi > mMaxLimit)
		{
			mMaxXi = mMaxLimit;
		}
	}
	else
	{
		mMinXi += deltaXi;
		if (mMinXi < -mMaxLimit)
		{
			mMinXi = -mMaxLimit;
		}
	}
}

bool Calculator::UpdateCache(float inX, float inY, float inVx, float inVy, 
	bool inCheckOccupied)
{
	mOverlapFlags = 0;
	mCachedVx = inVx;
	mCachedVy = inVy;
	
	int bigCenterX = int(kTwiceMagnification * inX);
	int bigCenterY = int(kTwiceMagnification * inY);

	int beginTwiceX, beginTwiceY, endTwiceX, endTwiceY, bigDx, bigDy;
	GetTwiceXYLimit(-inVy, inVx, bigCenterX, beginTwiceX, endTwiceX, bigDx);
	GetTwiceXYLimit(inVx, inVy, bigCenterY, beginTwiceY, endTwiceY, bigDy); 

	if (beginTwiceX >= endTwiceX || beginTwiceY >= endTwiceY)
	{
		mInfosEnd = mInfos;
		return false;
	}

	int rangeY = endTwiceY - beginTwiceY;
	mBufferKeeper1.NewBufferWithCapacity((endTwiceX - beginTwiceX) * rangeY, &mInfos);

	int bigVx = int(::floor(kMagnification * inVx + 0.5));
	int bigVy = int(::floor(kMagnification * inVy + 0.5));

	int xiStraight = int(- inVy * bigDx + inVx * bigDy);
	int eta = int(inVx * bigDx + inVy * bigDy);

	int deltaXi_y = bigVx;
	int deltaEta_y = bigVy;
	int deltaXi_x = -bigVy - bigVx * rangeY;
	int deltaEta_x = bigVx - bigVy * rangeY;

	float othersA = 0, othersB = 0;
	int othersCount = 0;
	int totalOthersXi = 0;
	float halfCurvature = 0.5f * mScaledCurvature;

	Info* p = mInfos;
	for (int x = beginTwiceX; x < endTwiceX; ++x)
	{
		for (int y = beginTwiceY; y < endTwiceY; ++y)
		{
			if (mMinEta <= eta && eta <= mMaxEta)
			{
				int xi = xiStraight - int((halfCurvature * eta) * eta);
				if (mMinXi <= xi && xi <= mMaxXi)
				{
					const SubpixelInfo* subPixelInfo = mSubpixelInfoProvider.GetInfo(x, y);
					if (0 == (subPixelInfo->occupiedState & (~kOccupiedBySelfMask)))
					{
						p->xi = xi;
						p->eta = eta;
						p->p = subPixelInfo;
						p->xiStraight = xiStraight;
						++p;
					}
					else if (inCheckOccupied)
					{
						char ch = subPixelInfo->occupiedState;
						mOverlapFlags |= ch;
						if (kOccupiedByOthersMask & ch)
						{
							++othersCount;
							totalOthersXi += xi;
							othersA += subPixelInfo->a;
							othersB += subPixelInfo->b;
						}
					}
				}
			}
			xiStraight += deltaXi_y;
			eta += deltaEta_y;
		}
		xiStraight += deltaXi_x;
		eta += deltaEta_x;
	}
	mInfosEnd = p;

	if (0 < othersCount &&
		mMinEta + othersCount * kMagnification >= mMaxEta)
	{
		int midGapXi = totalOthersXi / othersCount;
		int lowGap = MININT;
		int highGap = MAXINT;
		const Info* pp = mInfos;
		for (; pp < p; ++pp)
		{
			int xi = pp->xi;
			if (xi > midGapXi)
			{
				if (xi < highGap)
				{
					highGap = xi;
					if (lowGap + kMagnification >= highGap)
					{
						break;
					}
				}
			}
			else if (xi > lowGap)
			{
				lowGap = xi;
				if (lowGap + kMagnification >= highGap)
				{
					break;
				}
			}
		}
		if (pp == p && MININT != lowGap && MAXINT != highGap)
		{
			// if there is a gap among xi-coordinates,
			// it means that the current contour dot is placed
			// above the older stroke (which is quite narrow).
			return false;
		}
	}

	return mInfos + othersCount < p ||
		(othersA * (inVx * inVx - inVy * inVy) + 2.0f * othersB * inVx * inVy <= 0);
}

bool Calculator::CalculateF(float inA, float inB)
{
	mInfoWithMaxValue = NULL;
	float fmax = 0;

	for (Info* p = mInfos; p < mInfosEnd; ++p)
	{
		float vEta = 1 - mScaledCurvature * p->xiStraight;
		float vXi = mScaledCurvature * p->eta;
		float a = p->p->a;
		float b = p->p->b;

		float f = (vEta * vEta - vXi * vXi) * (a * inA + b * inB) +
			2 * vEta * vXi * (b * inA - a * inB);

		f /= (vEta *vEta + vXi* vXi);

		p->f = f;
		if (f > fmax)
		{
			fmax = f;
			mInfoWithMaxValue = p;
		}	
	}
	return (NULL != mInfoWithMaxValue);
}

void Calculator::SummarizeAB(float& outA, float& outB) const
{
	float a = 0, b = 0;
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		if (mMinXi <= p->xi && p->xi <= mMaxXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			a += p->p->a;
			b += p->p->b;
		}	
	}
	outA = a;
	outB = b;
}

static float FindPhi0AsRoot(float inTau)
{
	float phi = ::atan(::sqrt(inTau));
	float phi0 = phi;

	// limit the number of iterations
	int iterationReverseCounter = 10;
	while (--iterationReverseCounter)
	{
		float phi2 = 2 * phi;
		float s2 = ::sin(phi2);
		float c2 = ::cos(phi2);
		float s4 = (1 - c2) * (1 - c2);
		float F = 1.5f * inTau * (1 + c2) * (phi2 - s2) - phi * s4;
		float F1 = 6 * inTau * s2 * s2 - s4 - 
			phi2 * s2 * (3 * inTau + 2 - 2 * c2);
		float dPhi = F / F1;
		phi -= dPhi;
		if (::fabs(dPhi) < 0.0001f)
		{
			break;
		}
	}

	return (0 < iterationReverseCounter ? phi : phi0);
}

void Calculator::SummarizeABAndCalculateCurvature(float& outA, float& outB, 
		float& outCurvature)
{
	float sa, sb, momentA, momentB, count, sumOfEtaSquares;
	sa = sb = momentA = momentB = count = sumOfEtaSquares = 0.0f;
	
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		int eta = p->eta;
		if (mMinXi <= p->xi && p->xi <= mMaxXi &&
			mMinEta <= eta && eta <= mMaxEta)
		{
			float a = p->p->a;
			float b = p->p->b;
			float weight = a * a + b * b;

			sa += a;
			sb += b;
			momentA += eta * a;
			momentB += eta * b;
			sumOfEtaSquares += weight * float(eta) * eta;
			count += weight;
		}	
	}

	float curvature = 0;
	float q1 = sa * sa + sb * sb;
	if ( 0 < q1 && 0 < sumOfEtaSquares)
	{
		float p = (sa * momentB - sb * momentA) / q1;
		float r = p * count / sumOfEtaSquares;
		float t = r * p;

		if (kCurvatureCalculationFast == mCurvatureCalculationToken)
		{
			curvature = 0.5f * r / (1/11.0f +  ::sqrt(100/121.0f + 0.5f * t));
		}
		else if (kCurvatureCalculationBest == mCurvatureCalculationToken && 0 != t)
		{
			if (t < 0.1f)
			{
				curvature = 0.5f * r / (1 + 0.275f * t - 0.037276786f * t * t);
			}
			else
			{
				float phi0 = FindPhi0AsRoot(0.75f * t);
				float s = ::sin(phi0);
				float c = ::cos(phi0);

				curvature = r * 0.75f * c * (phi0 - s * c) / (phi0 * s * s);
			}
		}
	}
	mScaledCurvature = curvature;
	outCurvature = curvature * kTwiceMagnification;
	outA = sa;
	outB = sb;
}

void Calculator::SummarizeSideColors(float inXiShift, short int* outColorsComponents)
{
	float sumC[3], sumCF[3], sumF, sumFF;
	sumF = sumFF = 0.0f;
	for (int i = 0; i < 3; ++i)
	{
		sumC[i] = sumCF[i] = 0.0f;
	}

	int pixelsCountInt = 0;
	const float xiMax = ::sqrt(mSquaredLineWidth);
	const int safetyMargin = (max(mMaxXi, -mMinXi) + 
		max(mMaxEta, -mMaxEta)) / kMagnification;

	if (mSubpixelInfoProvider.IsSafeToGetUpperRightNeighbor(
			mInfos->p, safetyMargin))
	{
		for (const Info* p = mInfos; p < mInfosEnd; ++p)
		{
			if (mMinXi <= p->xi && p->xi <= mMaxXi &&
				mMinEta <= p->eta && p->eta <= mMaxEta)
			{
				float x = float(p->xi) / xiMax;
				if (::fabs(x) <= 1.0f)
				{
					++pixelsCountInt;
					float x2 = x * x;
					// "f" is an approximation of integral from function (1-x*x)^(3/2)
					float f = x * (1 - x2 * (0.491826706f -
						x2 * (0.045323591f + x2 * 0.035551738f)));
					sumF += f;
					sumFF += f * f;
					const float* colors = 
						mSubpixelInfoProvider.GetColorsAtGradientPoint(p->p);
					for (int i = 3; i--;)
					{
						// make the red color be stored the first
						// (in the SubpixelInfo::colors the order is: blue, green, red)
						const float& c = colors[2 - i];
						sumC[i] += c;
						sumCF[i] += c * f;
					}
				}
			}
		}
	}
	else
	{
		const float magnifiedXiShift = inXiShift * kTwiceMagnification;
		for (const Info* p = mInfos; p < mInfosEnd; ++p)
		{
			if (mMinXi <= p->xi && p->xi <= mMaxXi &&
				mMinEta <= p->eta && p->eta <= mMaxEta)
			{
				float x = (float(p->xi) - magnifiedXiShift) / xiMax;
				if (::fabs(x) > 1.0f)
				{
					continue;
				}

				++pixelsCountInt;
				// "f" is an approximation of integral from function (1-x*x)^(3/2)
				float x2 = x * x;
				float f = x * (1 - x2 * (0.491826706f -
						x2 * (0.045323591f + x2 * 0.035551738f)));
				sumF += f;
				sumFF += f * f;

				for (int i = 3; i--;)
				{
					// make the red color be stored the first
					// (in the SubpixelInfo::colors the order is: blue, green, red)
					const float& c = p->p->colors[2 - i];
					sumC[i] += c;
					sumCF[i] += c * f;
				}
			}	
		}
	}

	// calculate output, using statistics
	const float pixelsCount = float(pixelsCountInt);
	const float delta = sumFF * pixelsCount - sumF * sumF;
	for (int i = 0; i < 3; ++i)
	{
		if (0 != delta)
		{
			float c0 = sumC[i] * sumFF - sumCF[i] * sumF;
			float c1 = 0.58904862f * 
				(sumCF[i] * pixelsCount - sumC[i] * sumF);
			outColorsComponents[i] = (short int)((c0 - c1) / delta);
			outColorsComponents[i + 3] = (short int)((c0 + c1) / delta);
		}
		else
		{
			// unlikely to go here because there is great color gradient
			// at the place of contour dot.  As a precaution, show the average color
			outColorsComponents[i] = outColorsComponents[i + 3] = 
				(short int)(sumC[i] / pixelsCount);
		}
	}
}

void Calculator::AdjustXiBounds()
{
	if (NULL == mInfoWithMaxValue)
	{
		return;
	}
	int xiOfMax = mInfoWithMaxValue->xi;
	float threshold = mInfoWithMaxValue->f * mThresholdFactor;
	int lowXi = mMinXi - 1;
	int highXi = mMaxXi + 1;
	for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
	{
		int xi = p->xi;
		if (lowXi < xi && xi < highXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			if (p->f < threshold)
			{
				if (xi < xiOfMax)
				{
					lowXi = xi;
				}
				else
				{
					highXi = xi;
				}
			}
		}
	}
	mMinXi = lowXi + 1;
	mMaxXi = highXi - 1;
}

//////////////////////////////////////////////////////////////////////////////

static inline bool LessEqual(int inLeft, int inRight)
{
	return inLeft <= inRight;
}
static inline bool GreaterEqual(int inLeft, int inRight)
{
	return inLeft >= inRight;
}

inline Calculator::ConstInfoPtr* 
	Calculator::FilterAndSort(ConstInfoPtr* inBegin, ConstInfoPtr* inEnd, 
		int inXiLimit,
		CompareIntsFunction inCompareXi, CompareIntsFunction inCompareEta)
{
	ConstInfoPtr* qEnd = inBegin;
	for (const ConstInfoPtr* p = inBegin; p < inEnd; ++p)
	{
		int xi = (*p)->xi;
		if (inCompareXi(inXiLimit, xi))
		{
			continue;
		}

		int eta = (*p)->eta;
		ConstInfoPtr* q = inBegin;

		while (q < qEnd && inCompareEta(eta, (*q)->eta))
		{
			++q;
		}
		if (q > inBegin && inCompareXi(q[-1]->xi, xi))
		{
			continue;
		}
		if (q == qEnd)
		{
			*qEnd++ = *p;
			continue;
		}
		ConstInfoPtr* qSkipper = q;
		while (qSkipper < qEnd && inCompareXi(xi, (*qSkipper)->xi))
		{
			++qSkipper;
		}
		if (q > inBegin && eta == (*q)->eta)
		{
			--q;
		}
		// insert at q, then put values from range [qSkipper...qEnd[
		if (q == qSkipper)
		{
			for (ConstInfoPtr* r = qEnd++; r > q; --r)
			{
				r[0] = r[-1];
			}
			*q = *p;
		}
		else if (q + 1 == qSkipper)
		{
			*q = *p;
		}
		else
		{
			*q++ = *p;
			while (qSkipper < qEnd)
			{
				*q++ = *qSkipper++;
			}
			qEnd = q;
		}
	}
	return qEnd;
}


void Calculator::AdjustRectExceptFront()
{
	ConstInfoPtr *lowLimiters, *highLimiters;
	int onePartSize = mInfosEnd - mInfos;
	mBufferKeeper2.NewBufferWithCapacity(2 * onePartSize, &lowLimiters);
	highLimiters = lowLimiters + onePartSize;

	ConstInfoPtr *lowPtr = lowLimiters, *highPtr = highLimiters;

	const float threshold = mInfoWithMaxValue->f * mThresholdFactor;
	const int xiOfMax = mInfoWithMaxValue->xi;
	const int etaOfMax = mInfoWithMaxValue->eta;

	int xi_low = mMinXi - 1;
	int xi_high = mMaxXi + 1;
	
	// gather limiting points with eta > etaOfMax
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		if (p->f >= threshold)
		{
			continue;
		}

		int xi = p->xi; 
		if (p->eta < etaOfMax)
		{
			if (xi < xiOfMax) *lowPtr++ = p;
			else *highPtr++ = p;
		}
		else
		{
			if (xi < xiOfMax)
			{
				if (xi > xi_low) 
				{
					xi_low = xi;
				}
			}
			else if (xi < xi_high)
			{
				xi_high = xi;
			}
		}
	}
	// filter and sort lowLimiters
	ConstInfoPtr *lowEnd = FilterAndSort(lowLimiters, lowPtr, xi_low,
		GreaterEqual, LessEqual);
	ConstInfoPtr *highEnd = FilterAndSort(highLimiters, highPtr, xi_high,
		LessEqual, LessEqual);
	
	int savedXiLow = xi_low;
	int savedXiHigh = xi_high;
	int savedEtaLow = mMinEta - 1;
	long long maxArea = 0;

	bool savedIsEndToEnd = true;

	lowPtr = lowLimiters;
	highPtr = highLimiters;
	mIsEndToEnd = true;
	
	for (bool isEnd = false; !isEnd;)
	{
		int eta;
		if (lowPtr < lowEnd)
		{
			if (highPtr < highEnd && (*highPtr)->eta > (*lowPtr)->eta)
			{
				eta = (*highPtr)->eta;
				xi_high = (*highPtr)->xi;
				++highPtr;
			}
			else
			{
				eta = (*lowPtr)->eta;
				xi_low = (*lowPtr)->xi;
				++lowPtr;
			}
		}
		else
		{
			if (highPtr < highEnd)
			{
				eta = (*highPtr)->eta;
				xi_high = (*highPtr)->xi;
				++highPtr;
			}
			else
			{
				eta = mMinEta - 1;
				isEnd = true;
			}
		}

		long long area = (long long)(xi_high - xi_low) * (mMaxEta - eta);
		if (area > maxArea)
		{
			maxArea = area;
			savedXiLow = xi_low;
			savedXiHigh = xi_high;
			savedEtaLow = eta;
			savedIsEndToEnd = isEnd;
		}
	}
	mMinXi = savedXiLow + 1;
	mMaxXi = savedXiHigh - 1;
	mMinEta = savedEtaLow + 1;
	mIsEndToEnd = savedIsEndToEnd;
}

void Calculator::AdjustRectExceptBack()
{
	ConstInfoPtr *lowLimiters, *highLimiters;
	int onePartSize = mInfosEnd - mInfos;
	mBufferKeeper2.NewBufferWithCapacity(2 * onePartSize, &lowLimiters);
	highLimiters = lowLimiters + onePartSize;

	ConstInfoPtr *lowPtr = lowLimiters, *highPtr = highLimiters;

	const float threshold = mInfoWithMaxValue->f * mThresholdFactor;
	const int xiOfMax = mInfoWithMaxValue->xi;
	const int etaOfMax = mInfoWithMaxValue->eta;

	int xi_low = mMinXi - 1;
	int xi_high = mMaxXi + 1;
	
	// gather limiting points with eta > etaOfMax
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		if (p->f >= threshold)
		{
			continue;
		}

		int xi = p->xi; 
		if (p->eta > etaOfMax)
		{
			if (xi < xiOfMax) *lowPtr++ = p;
			else *highPtr++ = p;
		}
		else
		{
			if (xi < xiOfMax)
			{
				if (xi > xi_low) 
				{
					xi_low = xi;
				}
			}
			else if (xi < xi_high)
			{
				xi_high = xi;
			}
		}
	}
	// filter and sort lowLimiters
	ConstInfoPtr *lowEnd = FilterAndSort(lowLimiters, lowPtr, xi_low,
		GreaterEqual, GreaterEqual);
	ConstInfoPtr *highEnd = FilterAndSort(highLimiters, highPtr, xi_high,
		LessEqual, GreaterEqual);
	
	
	int savedXiLow = xi_low;
	int savedXiHigh = xi_high;
	int savedEtaHigh = mMaxEta + 1;
	long long maxArea = 0;

	int xi_low0 = xi_low;
	int xi_high0 = xi_high;

	bool savedIsEndToEnd = true;

	lowPtr = lowLimiters;
	highPtr = highLimiters;
	mIsEndToEnd = true;
	
	for (bool isEnd = false; !isEnd;)
	{
		int eta;
		if (lowPtr < lowEnd)
		{
			if (highPtr < highEnd && (*highPtr)->eta < (*lowPtr)->eta)
			{
				xi_high = (*highPtr)->xi;
				eta = (*highPtr)->eta;
				++highPtr;
			}
			else
			{
				xi_low = (*lowPtr)->xi;
				eta = (*lowPtr)->eta;
				++lowPtr;			
			}
		}
		else
		{
			if (highPtr < highEnd)
			{
				xi_high = (*highPtr)->xi;
				eta = (*highPtr)->eta;
				++highPtr;
			}
			else
			{
				eta = mMaxEta + 1;
				isEnd = true;
			}
		}
		long long area = (long long)(xi_high - xi_low) * (eta - mMinEta);
		if (area > maxArea)
		{
			maxArea = area;
			savedXiLow = xi_low;
			savedXiHigh = xi_high;
			savedEtaHigh = eta;
			savedIsEndToEnd = isEnd;
		}
	}
	mMinXi = savedXiLow + 1;
	mMaxXi = savedXiHigh - 1;
	mMaxEta = savedEtaHigh - 1;
	mIsEndToEnd = savedIsEndToEnd;
}

void Calculator::CutAtTheConcaveXi()
{
	if (NULL == mInfoWithMaxValue)
	{
		return;
	}
	int xiOfMax = mInfoWithMaxValue->xi;
	float fMax = mInfoWithMaxValue->f;

	// find the extreme subpixels with lowest and highest xi.
	int xiLeft = xiOfMax, xiRight = xiOfMax;
	float fLeft = fMax, fRight = fMax;
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		int xi = p->xi;
		if (mMinXi <= xi && xi <= mMaxXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			if (xi < xiLeft)
			{
				xiLeft = xi;
				fLeft = p->f;
			}
			else if (xi > xiRight)
			{
				xiRight = xi;
				fRight = p->f;
			}
		}
	}

	// find the points below the two lines f(xi) = fMax + k * (xi - xiOfMax);
	// that connect the point of maximum with the left and right points
	float leftK = xiLeft < xiOfMax ? (fMax - fLeft) / (xiOfMax - xiLeft) : 0;
	float rightK = xiRight > xiOfMax ? (fMax - fRight) / (xiOfMax - xiRight) : 0;
	int xiMin = mMinXi, xiMax = mMaxXi;

	// exclude some region around the subpixel with Max value
	int xiMin1 = (xiMin + xiOfMax) / 2;
	int xiMax1 = (xiMax + xiOfMax) / 2;

	float leftBottom = fMax, rightBottom = fMax;
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		int xi = p->xi;
		if (mMinXi <= xi && xi <= mMaxXi && 
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			if (xi < xiMin1)
			{
				float b = p->f + leftK * (xiOfMax - xi);
				if (b < leftBottom)
				{
					leftBottom = b;
					xiMin = xi;
				}
			}
			else if (xi > xiMax1)
			{
				float b = p->f + rightK * (xiOfMax - xi);
				if ( b < rightBottom)
				{
					rightBottom = b;
					xiMax = xi;
				}
			}
		}
	}
	mMinXi = xiMin;
	mMaxXi = xiMax;
}

bool Calculator::FindAlternativeXiBounds()
{
	if (NULL == mInfoWithMaxValue)
	{
		return false;
	}
	int xiOfMax = mInfoWithMaxValue->xi;

	ConstInfoPtr *pointers;
	mBufferKeeper2.NewBufferWithCapacity(mInfosEnd - mInfos, &pointers);
	int count = 0;
	if (0 < xiOfMax)
	{
		for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
		{
			int xi = p->xi;
			if (xi >= xiOfMax)
			{
				continue;
			}
			int a = 0, b = count;
			while (a < b)
			{
				int i = (a + b) >> 1;
				if (xi < pointers[i]->xi)
				{
					a = i + 1;
				}
				else 
				{
					b = i;
				}
			}
			if (a < count)
			{
				::memmove(&pointers[a+1], &pointers[a], (count - a) * sizeof(pointers[0]));
			}
			pointers[a] = p;
			++count;
		}
	}
	else // here xiOfMax < 0
	{
		for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
		{
			int xi = p->xi;
			if (xi <= xiOfMax)
			{
				continue;
			}
			int a = 0, b = count;
			while (a < b)
			{
				int i = (a + b) >> 1;
				if (xi > pointers[i]->xi)
				{
					a = i + 1;
				}
				else 
				{
					b = i;
				}
			}
			if (a < count)
			{
				::memmove(&pointers[a+1], &pointers[a], (count - a) * sizeof(pointers[0]));
			}
			pointers[a] = p;
			++count;
		}
	}

	float currentF = mInfoWithMaxValue->f + 1;
	int xi1OfLow = 0, xi2OfLow = 0, countLow = 0;
	int currentXi1 = 0, currentCount = 0;

	ConstInfoPtr currentMaxPtr = NULL;
	ConstInfoPtr maxPtr = NULL;

	float currentMaxF = 0;

	for (int i = 0; i < count; ++i)
	{
		ConstInfoPtr p = pointers[i];
		if (p->f < currentF)
		{
			if (currentCount > countLow)
			{
				countLow = currentCount;
				xi1OfLow = currentXi1;
				xi2OfLow = p->xi;
				maxPtr = currentMaxPtr;
			}
			currentMaxF = currentF = p->f;
			currentXi1 = p->xi;
			currentCount = 0;
		}
		else
		{
			++currentCount;
			if (p->f > currentMaxF)
			{
				currentMaxF = p->f;
				currentMaxPtr = p;
			}
		}
	}
	if (0 == countLow)
	{
		xi1OfLow = currentXi1;
		maxPtr = currentMaxPtr;
	}

	if (NULL != maxPtr && 3 < countLow + currentCount)
	{
		mInfoWithMaxValue = maxPtr;
		if (0 > xiOfMax)
		{
			mMinXi = xi1OfLow;
		}
		else
		{
			mMaxXi = xi1OfLow;
		}
	}

	return true;
}

bool Calculator::CheckAlignmentOfMaxPixel(float inReferencePointEta) const
{
	if (NULL == mInfoWithMaxValue)
	{
		return false;
	}

	float vEta = float(mInfoWithMaxValue->eta - inReferencePointEta * kTwiceMagnification);
	float vXi = float(mInfoWithMaxValue->xi);

	// compose a coefficient that is 1 for big distances but is smaller for
	// distances around the size of pixel (kTwiceMagnification).
	// This is necessary to suppress false calculations due to rounding (to pixel boundary) errors
	float absEta = ::fabs(vEta);
	float compensationCoef = absEta / (kTwiceMagnification + absEta);

	// rotate by the angle of curving
	float angle = 0.5f * mScaledCurvature * mInfoWithMaxValue->eta;
	vEta -= angle * vXi;
	vXi += angle * vEta;

	// this is the direction to the vertex of isosceles triangle 
	float vEta2 = vEta * vEta - vXi * vXi;
	float vXi2 = 2 * vXi * vEta;

	// rotate to XY
	float vx = vEta2 * mCachedVx - vXi2 * mCachedVy;
	float vy = vEta2 * mCachedVy + vXi2 * mCachedVx;

	float a = vx * vx - vy * vy;
	float b = 2 * vx * vy;
	float scalarProduct = a * mInfoWithMaxValue->p->a + b * mInfoWithMaxValue->p->b;
	float vectorProduct = b * mInfoWithMaxValue->p->a - a * mInfoWithMaxValue->p->b;

	return (compensationCoef * ::fabs(vectorProduct) < scalarProduct);
}

//-------------------------------------
bool Calculator::ApproximateByParabola(bool inCalculateAll)
{
	if (mRegardEtaVariations)
	{
		return ApproximateByParabola1(inCalculateAll);
	}

	float sx, sx2, sx3, sx4, sf, sxf, sx2f, sff, sy, syf;
	sx = sx2 = sx3 = sx4 = sf = sxf = sx2f = sff = sy = syf = 0;
	int s0 = 0;

	float k = kBestApproximationToParabolaFactor / mInfoWithMaxValue->f;

	for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
	{
		int xi = p->xi;
		if (mMinXi <= xi && xi <= mMaxXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			float f =  p->f;
			f -= k * f * f;

			++s0;
			sf += f;

			if (inCalculateAll)
			{
				sff += f * f;
				sy += p->eta;
				syf += f * p->eta;
			}

			sx += xi;
			sxf += (f *= xi);
			sx2f += f * xi;

			float xi2 = float(xi) * xi;
			sx2 += xi2;
			sx3 += xi2 * xi;
			sx4 += xi2 * xi2;
		}
	}
	if (s0 < 4)
	{
		return false;
	}

	float x22 = sx4 * s0 - sx2 * sx2;
	float x11 = sx2 * s0 - sx * sx;
	float x12 = sx3 * s0 - sx * sx2;
	float z1 = sxf * s0 - sx * sf;
	float z2 = sx2f * s0 - sx2 * sf;
	float Da = z2 * x11 - z1 * x12;
	float Db = z1 * x22 - z2 * x12;
	if (0 <= Da)
	{
		return false;
	}
	mDeltaXi = - 0.5f * Db / Da;
	if (::fabs(mDeltaXi) > float((mMaxXi - mMinXi) / 2))
	{
		return false;
	}

	if (!inCalculateAll)
	{
		// skip calculation of the other variables (besides mDeltaXi)
		return true;
	}
	
	if (NULL != mInfoWithMaxValue)
	{
		float dyf = sf * sy - syf * s0;
		mRegardEtaVariations = (::fabs(dyf) > 
			0.1f * s0 * (mMaxEta - mMinEta) * mInfoWithMaxValue->f);
	}

	float D = x11 * x22 - x12 * x12;
	float f11 = sff * s0 - sf * sf;
	if (0 >= D || 0 >= f11)
	{
		return false;
	}

	mParabolaA = Da / D;
	mParabolaB = Db / D;
	mParabolaC = (sf - mParabolaA * sx2 - mParabolaB * sx) / s0;

	mSquaredLineWidth = mDeltaXi * mDeltaXi - mParabolaC / mParabolaA;
	if (0 >= mSquaredLineWidth || mSquaredLineWidth >= mTooBigSquaredLineWidth)
	{
		return false;
	}

	mApproximationEpsilon = 1 + (mParabolaA * mParabolaA * x22 + 
		mParabolaB * mParabolaB * x11 + 
		2 * mParabolaA * mParabolaB * x12 -
		2 * mParabolaA * z2 -
		2 * mParabolaB * z1) / f11;

	return true;
}

///////////////////////////////////////////////////////////////////////////
float Calculator::GetLastApproximationEpsilon() const
{
	return mApproximationEpsilon;
}

void Calculator::ResetEtaCalculations()
{
	mRegardEtaVariations = false;
}
bool Calculator::IsEtaCalculations() const
{
	return mRegardEtaVariations;
}

///////////////////////////////////////////////////////////////////////////
OverrideToken Calculator::HandleOverrides(float inX0, float inY0)
{
	// check the subpixel that is located just under the maximum of parabola
	int twiceX = int(2 * inX0 - mCachedVy * mDeltaXi / kMagnification);
	int twiceY = int(2 * inY0 + mCachedVx * mDeltaXi / kMagnification);
	const SubpixelInfo* p = mSubpixelInfoProvider.GetInfo(twiceX, twiceY);

	const char occupiedState = p->occupiedState;
	if (0 == (occupiedState & kOutsideBoundsMask))
	{
		return (occupiedState & kOccupiedByOthersMask ? 
			kOverrideOtherStroke : kOverrideNone);
	}
	
	// adjust the parameters of the contour dot that goes beyond the picture bounds
	int xiOfMax = mInfoWithMaxValue->xi;
	float dxi = float(xiOfMax) - mDeltaXi;

	if (0 <= dxi)
	{
		// the "upper" part (above mDeltaXi) is knownint
		mDeltaXi += dxi * min(1.0f, dxi / mMaxLimit);
		mMinXi = int(xiOfMax - mDeltaXi) - 1;
		mMaxXi -= int(mDeltaXi);
	}
	else
	{
		// the "lower" part (below mDeltaXi) is known
		mDeltaXi += dxi * min(1.0f, -dxi / mMaxLimit);
		mMinXi -= int(mDeltaXi);
		mMaxXi = int(xiOfMax - mDeltaXi) + 1;
	}

	return kOverrideBeyondPictureBounds;
}

///////////////////////////////////////////////////////////////////////////
void Calculator::SetMaximumDotSize(float inSize)
{
	mMaxDotHalfsize = 0.5f * inSize;
	mTooBigSquaredLineWidth = 10 * mMaxDotHalfsize * mMaxDotHalfsize 
		* kTwiceMagnification * kTwiceMagnification;
}

void Calculator::SetCurvatureCalculationToken(int inToken)
{
	mCurvatureCalculationToken = inToken;
}

float Calculator::GetDeltaXi() const
{
	return mDeltaXi / kTwiceMagnification;
}

float Calculator::GetLineWidth() const
{
	return float(::sqrt(mSquaredLineWidth) * 0.375 * M_PI / kTwiceMagnification);
}

float Calculator::GetDotHalfsize() const
{
	float halfsize = ::sqrt(mSquaredLineWidth) / kTwiceMagnification;
	
	// clip the value to range (1..mMaxDotHalfsize)
	return max(1, min(halfsize, mMaxDotHalfsize));
}

bool Calculator::IsForwardOverlapping() const
{
	int midEta = (mMaxEta + mMinEta) / 2;
	for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
	{
		if (0 != (kOccupiedBySelfMask & p->p->occupiedState) &&
			midEta < p->eta && p->eta < mMaxEta &&
			mMinXi <= p->xi && p->xi <= mMaxXi)
		{
			return true;
		}
	}
	return false;
}

bool Calculator::IsBackwardOverlapping() const
{
	int midEta = (mMaxEta + mMinEta) / 2;
	for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
	{
		if (0 != (kOccupiedBySelfMask & p->p->occupiedState) &&
			mMinEta < p->eta && p->eta < midEta &&
			mMinXi <= p->xi && p->xi <= mMaxXi)
		{
			return true;
		}
	}
	return false;
}

void Calculator::NoteToAnchor()
{
	NoteToBooker(mAnchorDotBooker);
}

void Calculator::NoteToCurrent()
{
	NoteToBooker(mDotBooker);
}

void Calculator::MarkOccupiedTheCurrent()
{
	mCurrentDotBooker->AddUnmaskedReferencesTo(mTotalsBooker, 
		0, kOccupiedBySelfMask);
}

void Calculator::ChangeTheAnchorState(bool inOccupiedState)
{
	if (inOccupiedState)
	{
		mAnchorDotBooker.ChangeInfosOccupiedState(kOccupiedBySelfMask, 0);
	}
	else
	{
		mAnchorDotBooker.ChangeInfosOccupiedState(0, kOccupiedBySelfMask);
	}
}


void Calculator::UpdateBooking(char inSavingMask)
{
	mTotalsBooker.ChangeInfosOccupiedState(inSavingMask, kOccupiedBySelfMask);
	mTotalsBooker.RemoveAllReferences();
}

bool Calculator::IsMiddlePoint() const
{
	return mIsEndToEnd;
}

void Calculator::GetXiBounds(float& outMinXi, float& outMaxXi) const
{
	outMinXi = float(mMinXi) / kTwiceMagnification;
	outMaxXi = float(mMaxXi) / kTwiceMagnification;
}

float Calculator::GetMinEta() const
{
	return float(mMinEta) / kTwiceMagnification;
}

float Calculator::GetMaxEta() const
{
	return float(mMaxEta) / kTwiceMagnification;
}

//------------------------------------------------------------
bool Calculator::ApproximateByParabola1(bool inCalculateAll)
{
	float sx, sx2, sx3, sx4, sf, sxf, sx2f, sff, sy, syf, sxy, sx2y, sy2;
	sx = sx2 = sx3 = sx4 = sf = sxf = sx2f = sff = sy = syf = sxy = sx2y = sy2 = 0;
	int s0 = 0;

	float k = kBestApproximationToParabolaFactor / mInfoWithMaxValue->f;
	float theScaleToFitFloatRange = 1.0f / float(kTwiceMagnification);

	for (ConstInfoPtr p = mInfos; p < mInfosEnd; ++p)
	{
		int xi = p->xi;
		if (mMinXi <= xi && xi <= mMaxXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta)
		{
			float f =  p->f;
			f -= k * f * f;

			++s0;
			sf += f;

			if (inCalculateAll)
			{
				sff += f * f;
			}

			float eta = float(p->eta);
			float correctedXi = xi * theScaleToFitFloatRange;
			float xi2 = correctedXi * correctedXi;

			sy += eta;
			syf += f * eta;
			sxy += correctedXi * eta;
			sx2y += xi2 * eta;
			sy2 += eta * eta;

			sx += correctedXi;
			sxf += (f *= correctedXi);
			sx2f += f * correctedXi;

			
			sx2 += xi2;
			sx3 += xi2 * correctedXi;
			sx4 += xi2 * xi2;
		}
	}
	if (s0 < 4)
	{
		return false;
	}

	float x22 = sx4 * s0 - sx2 * sx2;
	float x11 = sx2 * s0 - sx * sx;
	float x12 = sx3 * s0 - sx * sx2;
	float z1 = sxf * s0 - sx * sf;
	float z2 = sx2f * s0 - sx2 * sf;
	float y1 = sxy * s0 - sx * sy;
	float y2 = sx2y * s0 - sx2 * sy;
	float q = sy2 * s0 - sy * sy;
	float p = syf * s0 - sy * sf;
	float Da_new = q * (x11 * z2 - x12 * z1) +
		p * (x12 * y1 - x11 * y2) + y1 * (z1 * y2 - z2 * y1);
	float Db_new = q * (z1 * x22 - z2 * x12) +
		p * (x12 * y2 - x22 * y1) + y2 * (z2 * y1 - z1 * y2);

	if (0 <= Da_new)
	{
		return false;
	}
	mDeltaXi = - 0.5f * Db_new / (Da_new * theScaleToFitFloatRange);
	
	if (::fabs(mDeltaXi) > float((mMaxXi - mMinXi) / 2))
	{
		return false;
	}

	if (!inCalculateAll)
	{
		// skip calculation of the other variables (besides mDeltaXi)
		return true;
	}

	float D = x11 * x22 - x12 * x12;

	float D_new = q * D + 2 * y1 * y2 * x12 -
		x11 * y2 * y2 - x22 * y1 * y1;

	
	float f11 = sff * s0 - sf * sf;
	if (0 >= D_new || 0 >= f11)
	{
		return false;
	}

	float Dd = p * D + y1 * (z2 * x12 - z1 * x22) +
		y2 * (z1 * x12 - z2 * x11);

	mParabolaA = Da_new / D_new;
	mParabolaB = Db_new / D_new;
	mParabolaD = Dd / D_new;
	mParabolaC = (sf - mParabolaA * sx2 - mParabolaB * sx - mParabolaD * sy) / s0;

	mSquaredLineWidth = mDeltaXi * mDeltaXi - 
		mParabolaC / (mParabolaA * theScaleToFitFloatRange * theScaleToFitFloatRange);
	if (0 >= mSquaredLineWidth || mSquaredLineWidth >= mTooBigSquaredLineWidth)
	{
		return false;
	}

	float u1 = mParabolaA * z2 + mParabolaB * z1 + mParabolaD * p;
	float u2 = mParabolaA * mParabolaA * x22 + mParabolaB * mParabolaB * x11 +
		mParabolaD * mParabolaD * q;
		
	float u3 = mParabolaA * mParabolaB * x12 + 
		mParabolaD * (mParabolaA * y2 + mParabolaB * y1);

	mApproximationEpsilon = 1 + (u2 + 2 * (u3 - u1)) / f11;
	mParabolaB *= theScaleToFitFloatRange;
	mParabolaA *= theScaleToFitFloatRange * theScaleToFitFloatRange;

	return true;
}

//-----------------------------------------------------------------
void Calculator::GetTwiceXYLimit(float inV1, float inV2, int inTwiceCoord,
	int& outMin, int& outMax, int& outFirstDelta)
{
	float xiF1 = mMinXi * inV1;
	float xiF2 = mMaxXi * inV1;

	float minValue = inTwiceCoord + min(xiF1, xiF2);
	float maxValue = inTwiceCoord + max(xiF1, xiF2);

	float etaF1 = mMinEta * inV2;
	float etaF2 = mMaxEta * inV2;

	if (float q = mScaledCurvature * inV1)
	{
		etaF1 += 0.5f * q * float(mMinEta) * mMinEta;
		etaF2 += 0.5f * q * float(mMaxEta) * mMaxEta;

		int etaOfExtremum = - int(inV2 / q);
		if (mMinEta < etaOfExtremum && etaOfExtremum < mMaxEta)
		{
			float valueAtExtremum = -0.5f * inV2 * inV2 / q;
			if (0 < q)
			{
				// it is miminum at the extremum
				minValue += valueAtExtremum;
				maxValue += max(etaF1, etaF2);
			}
			else
			{
				// it is maximum at the extremum
				minValue += min(etaF1, etaF2);
				maxValue += valueAtExtremum;
			}
		}
		else
		{
			minValue += min(etaF1, etaF2);
			maxValue += max(etaF1, etaF2);
		}

	}
	else
	{
		minValue += min(etaF1, etaF2);
		maxValue += max(etaF1, etaF2);
	}

	outMin = int(::floor(minValue / kMagnification));
	outMax = int(::floor(maxValue / kMagnification)) + 1;
	outFirstDelta = (outMin * kMagnification) - inTwiceCoord;
}

void Calculator::NoteToBooker(Booker& ioBooker)
{
	mCurrentDotBooker = &ioBooker;
	ioBooker.RemoveAllReferences();
	for (const Info* p = mInfos; p < mInfosEnd; ++p)
	{
		if (mMinXi <= p->xi && p->xi <= mMaxXi &&
			mMinEta <= p->eta && p->eta <= mMaxEta &&
			0 == ((kOccupiedBySelfMask | kOutsideBoundsMask) & p->p->occupiedState))
		{
			ioBooker.AddReference(const_cast<SubpixelInfo*>(p->p));
		}	
	}	
}